#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>

#include "Vector/vector.03.h"
#include "Color/color.01.h"
#include "Timer/timer.01.h"
#include <vector>
#include <deque> 

#define SQR(x) ((x)*(x))

#include "ics.h"

using namespace std;

typedef unsigned int uint;
typedef long unsigned int luint;
typedef double real;



const double kB = 1.3806488e-23;

const int gsizex = 80;
const int gsizey = 80;
const int gsizez = 500;

double gstepx = 3;
double &gstepy = gstepx;
double &gstepz = gstepx;



double drij2max = SQR(3.0); 
double dt = 1e-3;
uint period = 15;
uint nsteps = period/dt;

double lU = 1e-9; // m [1nm]
double eU = 1.60217646e-19; // J [1ev]
double mU = 1.660538921e-27; // kg [amu]
double tU = lU*sqrt(mU/eU); // sec
double vU = lU/tU; // sec

double g = -9.81; // sec

luint memmax = 500; //[MB]	

bool ssF = false;
double ss_period = 0.01;

bool coolingF = false;

double all_memsize;

int ii(1), e00(10000), e11(200), e01(2000);


struct P {
	Vector<real> R;
	Vector<real> V;
	Vector<real> F;
	int type;
};
vector<P> p;
vector<P*> grid[2*gsizex][2*gsizey][2*gsizez];

struct PW {
	double e;
	double s; 
	double rc;
};
PW pw[10][10];

struct PT {		
	double m;
	Color color;	
};
vector< PT > pt;		

uint ntot;
				
Timer tr, trf;
		
double rnd() {
	return (double)rand()/RAND_MAX;
};

int p_size;

template <int ROW, int COL>
void packToHCP(vector<P> &p, const PW (&pw)[ROW][COL], int type,  Vector<double> R, Vector<int> S, Vector<double> V, double t) {
			
	double dr = pw[type][type].s*0.977;
			
	for (int i=0; i<S.x; i++)
		for (int j=0; j<S.y; j++) 
			for (int k=0; k<S.z; k++){
					
				P pi;	
					
				pi.R.x = dr*i + (j%2==0 ? 0: 0.5*dr) + (k%2==0 ? 0: 0.5*dr);
				pi.R.y = 0.5*sqrt(3)*dr*j + (k%2==0 ? 0: sqrt(3.0)/6.0*dr);
				pi.R.z = sqrt(6.0)/3.0*dr*k;
					
				pi.V = Vector<double>(t*(1-2*rnd()), t*(1-2*rnd()), t*(1-2*rnd()));

				pi.R  += R;
						
				pi.V  += V;
				pi.type = type;	
				p.push_back(pi);	
					
			}
}


template <int ROW, int COL>
void packToHCP(vector<P> &p, const PW (&pw)[ROW][COL], int type1, int type2, Vector<double> R, Vector<int> S, Vector<double> V, double t) {
			
	double dr = pw[type1][type1].s*0.977;
			
	for (int i=0; i<S.x; i++)
		for (int j=0; j<S.y; j++) 
			for (int k=0; k<S.z; k++){
					
				P pi;	
					
				pi.R.x = dr*i + (j%2==0 ? 0: 0.5*dr) + (k%2==0 ? 0: 0.5*dr);
				pi.R.y = 0.5*sqrt(3)*dr*j + (k%2==0 ? 0: sqrt(3.0)/6.0*dr);
				pi.R.z = sqrt(6.0)/3.0*dr*k;
					
				pi.V = Vector<double>(t*(1-2*rnd()), t*(1-2*rnd()), t*(1-2*rnd()));
				
				pi.R.x-= 0.5*S.x*dr;
				pi.R.y-= 0.5*S.y*0.5*sqrt(3)*dr;
				pi.R.z-= 0.5*S.z*sqrt(6.0)/3.0*dr;
				
				pi.R  += R;
						
				pi.V  += V;
				
				pi.type =  (i==0 || j==0 || k==0 || i==S.x-1 || j==S.y-1 || k==S.z-1) ? type2 : type1;
				
				p.push_back(pi);	
					
			}
}


inline void calc_force(vector<P*> &pi) {
	
	
	for (int i=0;i<pi.size()-1;i++)
		for (int j=i+1;j<pi.size();j++) {
						
			real dx = pi[i]->R.x - pi[j]->R.x; // 1
			real dy = pi[i]->R.y - pi[j]->R.y; // 1
			real dz = pi[i]->R.z - pi[j]->R.z; // 1
								
			real r2 = dx*dx + dy*dy + dz*dz; // 5

			if ( r2 > pw[pi[i]->type][pi[j]->type].rc ) continue; 
							
			real r4 	= r2*r2; // 1
			real r6i 	= 1.0/(r2*r4); // 2
			real r12i = r6i*r6i;  // 1
			
			real f = pw[pi[i]->type][pi[j]->type].e*r6i*r4*(r12i-r6i); // 4
							
			real dxf = dx*f;
			real dyf = dy*f;
			real dzf = dz*f;
							
			pi[i]->F.x += dxf; // 2
			pi[i]->F.y += dyf; // 2
			pi[i]->F.z += dzf; // 2
						
			pi[j]->F.x -= dxf; // 2
			pi[j]->F.y -= dyf; // 2
			pi[j]->F.z -= dzf; // 2
							
	}
}


inline void calc_force(vector<P*> &pi,  vector<P*> &pj) {
					
	for (int i=0;i<pi.size();i++)
		for (int j=0;j<pj.size();j++) {
						
						
			real dx = pi[i]->R.x - pj[j]->R.x; // 1
			real dy = pi[i]->R.y - pj[j]->R.y; // 1
			real dz = pi[i]->R.z - pj[j]->R.z; // 1
								
			real r2 = dx*dx + dy*dy + dz*dz; // 5

			if ( r2 > pw[pi[i]->type][pj[j]->type].rc ) continue; 
							
			real r4 	= r2*r2; // 1
			real r6i  = 1.0/(r2*r4); // 2
			real r12i = r6i*r6i;  // 1
			
			real f = pw[pi[i]->type][pj[j]->type].e*r6i*r4*(r12i-r6i); // 4
			
			real dxf = dx*f;
			real dyf = dy*f;
			real dzf = dz*f;
			
			pi[i]->F.x += dxf; // 2
			pi[i]->F.y += dyf; // 2
			pi[i]->F.z += dzf; // 2
			
			pj[j]->F.x -= dxf; // 2
			pj[j]->F.y -= dyf; // 2
			pj[j]->F.z -= dzf; // 2
		}
	
					
}

double min(double value1, double value2, double value3, double value4) {
	
	return min(min(min(value1, value2), value3), value4);
	
}


void save_ss(double t) {
	
	char foldername[256];
	sprintf(foldername, "%i.%i.%i.%i", ii, e00, e11, e01);
	mkdir(foldername, 0777);
			
	char filename[256];
	sprintf(filename, "%s/ss.%.3f.dat", foldername, t);
		
	fprintf(stderr,"%s\n", filename);
			
	FILE *fp = fopen (filename , "wb");
		
	struct MEM { Vector<float> R; float type;} *mem;
	mem = new MEM[p_size];
	
	for (int i=0; i<p_size; i++) {
		mem[i].R = p[i].R;
		mem[i].type = p[i].type;
	}
		
	fwrite(mem, sizeof(MEM), p_size, fp);
	fclose(fp);
	delete [] mem;
	
}


void run(){
					
	for(int i=0; i<p_size; i++) p[i].F*=0;

	double t(0);
	
	
	
	double glenx = gsizex*gstepx;
	double gleny = gsizey*gstepy;
	double glenz = gsizez*gstepz;
	
	double g2sizex = 2*gsizex;
	double g2sizey = 2*gsizey;
	
	
	double dt = 0.00001/tU; // in 2 powers ?? ?
	
	double cdt = 1.1;
	double dtmax = 1e-3/tU;
	
	double dRmax = 1.0/100;
	double dEmax = pw[1][1].e*5; // minimum potential well depth 
	
	
	double Vmax(0);
	double Fmax(0);
	
	int step(0);
	
	save_ss(t);
	double ss_t = ss_period = 0.001/tU;
	
	do {
		
		
		double dt05 = 0.5*dt;
		
		tr.set();
		//printf("%i: #0\n", t);				
		for (int xi=0; xi<g2sizex; xi++)
		for (int yi=0; yi<g2sizey; yi++)
		for (int zi=0; zi<gsizez; zi++)
			grid[xi][yi][zi].resize(0);
		
		//piano_hilbert_ordering();
		
		//printf("%i: #1\n", t);				
		
		for(int i=0; i<p_size; i++) {
			
			double dt05_over_m = dt05/pt[p[i].type].m;
			p[i].R+= (p[i].V +  dt05_over_m*p[i].F)*dt;
			p[i].V+= dt05_over_m*p[i].F;
		
			// boundaries	
			if (p[i].R.z < 0) {
				p[i].R.z*=-1;  p[i].V.z*=-1; 
			} 
			
			if (p[i].R.x < -glenx || +glenx < p[i].R.x || p[i].R.y < -gleny || +gleny < p[i].R.y || glenz < p[i].R.z) {
				//printf("%i:( %g %g %g) out of range (%i %i)\n", i, p[i].R.x, p[i].R.y, p[i].R.z, p_size,  p.size());				
				p[i--] = p[--p_size];
				
				continue;
			} 
		
			// greed update
			grid[int(p[i].R.x/gstepx)+gsizex][int(p[i].R.y/gstepy)+gsizey][int(p[i].R.z/gstepz)].push_back(&p[i]);
			
			// forces
			p[i].F.x=p[i].F.y=p[i].F.z=0;
		}
						
		trf.set();

		//printf("%i: #2\n", t);				
		for (int xi=0; xi<g2sizex; xi++)
		for (int yi=0; yi<g2sizey; yi++)
		for (int zi=0; zi<gsizez; zi++) {
							
			vector<P*> &pi = grid[xi][yi][zi];
			
			if (pi.size()==0) continue;
							
			calc_force(pi);
			
			if (0 < xi && 0 < yi && zi+1<gsizez)					calc_force(pi, grid[xi-1][yi-1][zi+1]); // 12
			if (0 < xi && zi+1<gsizez)									calc_force(pi, grid[xi-1][yi  ][zi+1]); // 1			
			if (0 < xi && yi+1<g2sizey)								calc_force(pi, grid[xi-1][yi+1][zi  ]); // 13
			if (0 < xi && yi+1<g2sizey && zi+1<gsizez)			calc_force(pi, grid[xi-1][yi+1][zi+1]); // 6
			
			if (0 < yi && zi+1<gsizez)									calc_force(pi, grid[xi  ][yi-1][zi+1]); // 4
			if (zi+1<gsizez)												calc_force(pi, grid[xi  ][yi  ][zi+1]); // 2
			if (yi+1<g2sizey && zi+1<gsizez)							calc_force(pi, grid[xi  ][yi+1][zi+1]); // 5
			if (yi+1<g2sizey)												calc_force(pi, grid[xi  ][yi+1][zi  ]); // 8
			
			if (xi+1<g2sizex && 0 < yi && zi+1<gsizez)			calc_force(pi, grid[xi+1][yi-1][zi+1]); // 11
			if (xi+1<g2sizex && zi+1<gsizez)							calc_force(pi, grid[xi+1][yi  ][zi+1]); // 3
			if (xi+1<g2sizex )											calc_force(pi, grid[xi+1][yi  ][zi  ]); // 10
			if (xi+1<g2sizex && yi+1<g2sizey && zi+1<gsizez)	calc_force(pi, grid[xi+1][yi+1][zi+1]); // 7
			if (xi+1<g2sizex && yi+1<g2sizey)						calc_force(pi, grid[xi+1][yi+1][zi  ]); // 9
							
		}
		trf.get();

		
		double V2max(0), F2max(0);
		for(int i=0; i<p_size; i++) {
			
			p[i].F.z += g*pt[p[i].type].m;
			p[i].V+= p[i].F*(dt05/pt[p[i].type].m);
			
			if (p[i].V.sqr() > V2max) V2max = p[i].V.sqr();
			if (p[i].F.sqr() > F2max) F2max = p[i].F.sqr();
			
		}
		
		t+=dt;
		
		tr.get();		
		
		if ( fabs(t - ss_t) < 1e-12 ) { 
			save_ss(t*tU);
			ss_t+=ss_period;
			dt = dtmax;
		}
		
		step++;
		
		Vmax = sqrt(V2max);
		Fmax = sqrt(F2max);
		
		dt = min(dRmax/Vmax, dEmax/(Fmax*Vmax), cdt*dt, dtmax);

		if ( ss_t < t+dt) dt = ss_t-t;
 		
		
		//fprintf(stderr,"%9i  %10.5e  %10.5e  %10.2f  %10.2e %10.2e",  step, t,  dt, Vmax, p_size/tr.sin, tr.sin_tot);
		fprintf(stdout,"%9i  %10.5e    %8.2e = min ( %8.2e || %8.2e )    %10.2e, %10.2f",  step, t,  dt, dRmax/Vmax, dEmax/(Fmax*Vmax), p_size/tr.sin, tr.sin_tot);
		fprintf(stdout,"\n");
		//fprintf(stderr,"\n");
		fflush(stdout);				
		//fprintf(stderr,"%9i %10.4f %6.4f %9i %10.4f %10.4f %4.1f%%   %10.2e  %10.4f  %10.4f\n",t, t*tU, dt, p_size, trf.sin, tr.sin,  100*(tr.sin-trf.sin)/tr.sin,  p_size/tr.sin, 1/tr.sin, all_memsize);
						
	} while (t<period);
	
}


int main(int argc, char *argv[]) {

	
	ICs ics;
	ics.make_WTC();
	
	
	e00 = 250000;
	e11 =    200;
	e01 =  20000;
	
	ics.pw[0][0].e *= e00; 
	ics.pw[1][1].e *= e11;
	ics.pw[0][1].e = ics.pw[1][0].e *= e01;
	
	p.resize(ics.p_size);
	p_size = ics.p_size;
	for (int i=0; i<ics.p_size; i++) {
		p[i].R = ics.p[i].R;
		p[i].V = ics.p[i].V;
		p[i].type = ics.p[i].type;
	}
	ics.p.resize(0);

	// copy PW data 
	for (int i=0; i<10; i++)
	for (int j=0; j<10; j++) {
		pw[i][j].e = ics.pw[i][j].e;
		pw[i][j].s = ics.pw[i][j].s;
		pw[i][j].rc = ics.pw[i][j].rc;
	}
	
	// copy PT data 
	pt.resize(ics.pt.size());
	for (int i=0; i<pt.size(); i++)
		pt[i].m = ics.pt[i].m;
	
	eU = ics.eU;
	vU = ics.vU;
	tU = ics.tU;
	g = -9.81*SQR(ics.tU)/ics.lU; 	

	
	period = 4.0/tU; //sec
	gstepx = 3;
	ss_period = 0.01/tU;  //sec
	
	run();

	return 0;
	
}
